/*Group 04
José Renato da Silva Bonatelli 7288964
Matheus Martins Teixeira 7277482
Rafael de Paula Pinto Paravia 71522567*/

/*Simulating the Black Scholes Model with a parallel approach of the Monte Carlo algorithm*/

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <pthread.h>
#define MAX(A,B) A>B?A:B
#define NUM_THREADS 50

double s,e,r,v,t,mean;
double *trials;
double threadSum[NUM_THREADS];
double threadSum2[NUM_THREADS];
unsigned int M;

/*Thread function*/
void *blackScholes(void *i){
        unsigned int k = (unsigned int) i;
        /*Each thread calculates a fraction of the total value*/
        unsigned int begin = k*M/NUM_THREADS;
        unsigned int end = (k==(NUM_THREADS-1))?M:(k+1)*M/NUM_THREADS;
        unsigned int j;
        double temp,sum;
        const unsigned int A = 1664525, C = 1013904223;
        unsigned int rr= time(NULL) +k; /*Initial seed*/
        double aux = (r-v*v*0.5)*t;
        double aux2 = v*sqrt(t);
        double aux3 = exp(-r*t);
        sum  = 0;
        for (j = begin; j < end; j++){
                rr = rr*A+C;        /*Calculating a pseudo-random value*/
                temp = s*exp(aux+aux2*((float)rr/UINT_MAX));
                trials[j] = aux3 * (MAX(temp-e,0));
                sum  += trials[j]; /*Partial sum*/

        }
        threadSum[k] = sum;         /*Storing the partial sum*/
        pthread_exit(NULL);
}

/*Function to calculate the standard deviation*/
void *stdDevCalc(void *i){
        unsigned int k = (unsigned int) i;
        unsigned int begin = k*M/NUM_THREADS;
        unsigned int end = (k==(NUM_THREADS-1))?M:(k+1)*M/NUM_THREADS;
        unsigned  int j;
        double sum=0;
        double mean2 = mean;
        for (j = begin; j < end; j++)
                sum  += pow(trials[j]-mean2,2);
        threadSum2[k] = sum;
        pthread_exit(NULL);
}

int main(void){
        int i, returnValue;
        double total=0;
        pthread_t threads[NUM_THREADS];
        scanf(" %lf %lf %lf %lf %lf %d ",&s,&e,&r,&v,&t,&M);
        /*Input:
            s = Stock Price
            e = Strike Price
            r = Interest rate
            v = Volatility
            t = Time to expect
            m = Number of iterations */
        trials = (double*) malloc(M*sizeof(double));
        if(trials == NULL){printf("Error on memory allocation\n");exit(-1);} /*Error checking*/
        for(i = 0; i < NUM_THREADS; i++){
                returnValue = pthread_create(&threads[i], NULL, blackScholes, (void *)i);/*Starting the threads*/
                if(returnValue){    /*Error checking*/
                        printf("Error on pthread_create : error number %d\n", returnValue);
                        exit(-1);
                }
        }
        for(i = 0; i < NUM_THREADS; i++){
                pthread_join(threads[i], NULL); /*Waiting for all threads to finish*/
                total += threadSum[i];
        }
        mean = (double)total/M;             /*Calculating the mean value*/
        total = 0;
        for(i = 0; i < NUM_THREADS; i++){   /*Using threads to calculate the standard deviation*/
                returnValue = pthread_create(&threads[i], NULL, stdDevCalc, (void *)i);
                if(returnValue){
                        printf("Error on pthread_create : error number %d\n", returnValue);
                        exit(-1);
                }
        }
        for(i = 0; i < NUM_THREADS; i++){
                pthread_join(threads[i], NULL); /*Waiting for all threads to finish*/
                total += threadSum2[i];
        }
        double stddev = sqrt((double)total/M); /*Calculating the standard deviation*/

        double confwidth = 1.96 * stddev / sqrt(M); /*Confidence interval*/
        double confmin = mean - confwidth;          /*Maximum value*/
        double confmax = mean + confwidth;          /*Minimum value*/
        /*Output*/
        printf("S        %lf\nE        %lf\nr        %lf\nsigma    %lf\nT        %lf\nM        %d\n",s,e,r,v,t,M);
        printf("Confidence interval :  (%lf,%lf)\n",confmin,confmax);
        pthread_exit(NULL);
        free(trials);
        return 0;

}

